FUNCTION calc_stock_recruitment
  {
	/*
	The following code is used to derive unfished
	spawning stock biomass bo and the stock-
	recruitment parameters for the:
	Beverton-Holt Model
		-Rt=k*Ro*St/(Bo+(k-1)*St)*exp(delta-0.5*tau*tau)
	Ricker Model
		-Rt=so*St*exp(-beta*St)*exp(delta-0.5*tau*tau)
		
	Dec 6, 2010.  Modified code to allow empirical weight-at-age data
				This required a fecundity-at-age matrix.  Need to 
				project spawning biomass into the future.
	CHANGED bo should be a function of the average natural mortality
	TODO update phib calculation if age-specific M's are used.
	
	May 6, 2010.  Changed phib calculation based on average M 
	in cases where M varies over time. Otherwise bo is biased based
	on the initial value of M.
	
	CHANGED Need to adjust spawning biomass to post fishery numbers.
	CHANGED Need to adjust spawners per recruit (phib) to average fecundity.
	*/ 
	int i,itmp;
	dvariable tau = (1.-prho)/varphi;
	dvar_vector tmp_rt(2,nyr);
	dvar_vector lx(1,s_nages); lx=1;
	for(i=1+1;i<=s_nages;i++) lx(i)=lx(i-1)*exp(-m);
	//lx(nage)/=(1.-exp(-m_bar));
	
	//dvariable phib = (lx*exp(-m*cntrl(13))) * avg_fec;	//SM Dec 6, 2010
	dvariable phib = (lx*exp(-m)) * s_w_age*s_mat_age;	//SM Dec 6, 2010

	dvariable so = alpha/phib;		//max recruits per spawner
	dvariable beta;
	bo = ro*phib;  					//unfished spawning biomass
	
	//sbt=rowsum(elem_prod(N,fec));			//SM Dec 6, 2010
	//CHANGED adjusted spawning biomass downward by ctrl(12)
	//AQUI HAY QUE CAMBIAR Y TRAER LA BIOM DESOV POR AÑO desde sbiom y usando peak spaw
	//OJO REVISAR tspawn y que hay delay en un year
	for(i=1;i<=nyr;i++)
	{
		itmp = tspawn(i);
		ssb(i)=sbiom(itmp);
	}

	//for(i=syr;i<=nyr;i++)
	//{
	//	sbt(i) = elem_prod(N(i),exp(-Z(i)*cntrl(13)))*fec(i);
	//}
	//sbt(nyr+1) = N(nyr+1)*fec(nyr+1);
	//cout<<"sbt\n"<<sbt<<endl;
	//exit(1);
	//dvar_vector tmp_st=sbt(syr,nyr-sage).shift(syr+sage);
	
	dvar_vector tmp_st(1,nyr-1);
	for(i=1;i<nyr;i++)
	{
		tmp_st(i) = ssb(i);
	}
	
	switch(int(cntrl(2)))
	{
		case 1:
			//Beverton-Holt model
			beta = (alpha-1.)/bo;
			for(i=1;i<nyr;i++)
			{
				//tmp_rt=elem_div(so*tmp_st,1.+beta*tmp_st);
				tmp_rt(i+1)=so*tmp_st(i)/(1.+beta*tmp_st(i));
			}
			break;
		case 2:
			//Ricker model
			//dvariable so=alpha/phib;
			beta=log(alpha)/bo;
			for(i=1;i<nyr;i++)
			{
				//tmp_rt=elem_prod(so*tmp_st,exp(-beta*tmp_st));
				tmp_rt(i+1)=so*tmp_st(i)*mfexp(-beta*tmp_st(i)));
			}

		break;
	}
	
	//residuals in stock-recruitment curve
	rt=exp(log_rt(1+s_nages,nyr)); //trans(N)(1)(syr+1,nyr);
	delta = log(rt)-log(tmp_rt)+0.5*tau*tau; //residuals from S-R
	
	if(verbose)cout<<"**** Ok after calc_stock_recruitment ****"<<endl;
	
  }
